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Abstract 

The following question is addressed: under what conditions can 
a strange diffusive process, defined by a semi-dynamical V-Langevin 
equation or its associated Hybrid kinetic equation (HKE), be described 
by an equivalent purely stochastic process, defined by a Continuous 
Time Random Walk (CTRW) or by a Fractional Differential Equa- 
tion (FDE)? More specifically, does there exist a class of V-Langevin 
equations with long-range (algebraic) velocity temporal correlation, 
that leads to a time-fractional superdiffusive process? The answer 
is always affirmative in one dimension. It is always negative in two 
dimensions: any algebraically decaying temporal velocity correlation 
(with a Gaussian spatial correlation) produces a normal diffusive pro- 
cess. General conditions relating the diffusive nature of the process 
to the temporal exponent of the Lagrangian velocity correlation (in 
Corrsin approximation) are derived. 



1 Introduction 

Strange Transport has been the object of intense studies in recent years 
(A very recent qualitative review is [1]). (We use the terminology "Strange 
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transport" [2] rather than "Anomalous transport" which is customary in 
Dynamical Systems theory, but has a different meaning in Plasma transport 
theory). The concept first appeared in the theory of stochastic processes, 
especially the theory of Continuous Time Random Walks (CTRW) [3], [4]. 

Consider a disordered system (e.g. a turbulent fluid or plasma). The 
position x(t) at time t of one of its particles is determined by its interac- 
tions with the other particles and/or with external sources. In a strongly 
disordered system these interactions are modelled by a random field. The 
statistical description of the system thus involves the definition of an ensem- 
ble of realizations. One of the important quantities describing the system 
is the mean square deviation (MSD) of the random variable x(t). In many 
cases, this is represented asymptotically (t — > oo) by a simple increasing 
function of time: 

(5x 2 {t)) = M t", (1) 

where 5x(i) = x(i) — (x(t)) . The value of the "diffusion exponent" ji 
characterizes the diffusion regime of the system 1 : 



Diffusion exponent 


Regime 


< /x < 1 


SUBDIFFUSIVE 


fi = l 


(NORMAL) DIFFUSIVE 


1 < /! < 2 


SUPERDIFFUSIVE 



(2) 



The subdiffusive and the superdiffusive regimes are the ones called " STRANGE" . 
It will be seen in the forthcoming sections that the concept of " strange trans- 
port" involves much more than a simple statement about the behaviour of 
the MSD. 

All present theories of strange transport are of stochastic nature. From the 
very voluminous literature we only cite here a few of the more comprehensive, 
physically oriented review papers and books: [3] - [13], where additional 
references will be found. 

As stated in [2], transport theories (e.g., for fluids or plasmas) can be 
constructed on three levels. 

a) A purely statistical mechanical theory would be based on the kinetic 
equation for a set of interacting charged particles, combined with Maxwell's 
equations. This would be the most fundamental description, but becomes 
impossibly complicated in practice. 

f3) The next level of description would be a compromise, based on a "semi- 
dynamical" model of particles moving according to the laws of mechanics 

1 There exist superdiffusive regimes with \i > 2, but they wiil not be considered here. 
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(Newton's equation), but under the action of a fluctuating field, representing 
the action of the turbulent environment. This leads to stochastic ordinary 
differential equations of the Langevin type. More specifically, we consider 
V-Langevin equations [2], i.e. equations for the position of a "particle": 

^ = v[x(f),f], x(0)=0 (3) 

where the right hand side represents a fluctuating velocity field, defined 
by its statistical properties. We always consider a divergenceless velocity 
field: 



V-v(x,t) = 0. (4) 

Associated with (3) is a kinetic-type first-order partial differential equa- 
tion for the fluctuating particle distribution function /(x, t) (whose average 
is the density profile), called a hybrid kinetic equation (HKE) [2]. By 
definition, the characteristics of this equation are precisely the V-Langevin 
equations (3): 

a t /(x,*)+v(x,*)-V/(x,*) = 0. (5) 
Due to the property (4) this equation is equivalent to the following: 

t /(x,t) + V-[v(x,t) /(x,i)] = 0. (6) 

This equation ensures the conservation of the number of particles, and is 
rightly interpreted as a kinetic equation. 

7) The last level of description is the level of continuous time random 
walk (CTRW) theories and fractional diffusion equations (FDE), in 

which the deterministic dynamical laws are completely given up, and replaced 
by a purely random process. 

All these concepts will be defined and used in the forthcoming sections. 
These three levels are, of course, interrelated: each one of them should be 
justifiable as an approximation of the more fundamental one. 

The strange transport theories associated with CTRW and with FDE 
have been very successful in recent years in modelling many peculiar behav- 
iors observed, among other fields, in plasma and fusion physics. There is, 
however, still a serious gap in justifying these purely stochastic models on 
the basis of a molecular theory, i.e., theories of the type a), or at least (5). 
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In [2] we treated several models exhibiting strange transport, showing that 
a stochastic CTRW can be associated with a semi-dynamical V-Langevin 
equation under certain limiting conditions. All these models (diffusion in a 
fluctuating electrostatic field or in a fluctuating magnetic field) were either 
diffusive, or subdiffusive. One of the purposes of the present work is to try to 
determine whether there exist V-Langevin equations leading to temporal su- 
perdiffusive behavior and associating them with a CTRW or with a fractional 
diffusion equation. 

A first work following the same philosophy was done by West et al. [14]. 
In that paper a specially simplified model of a one-dimensional, two-state 
CTRW was used. This means that a "particle" performs a CTRW in one 
dimension, moving with a velocity that is constant in absolute value: ±V. 
The particle moves during a given (random) time, after which it suddenly 
jumps and reverses its velocity. The velocity autocorrelation function is given 
as an algebraically decaying function of time. In a first step of their work, 
the jumps are supposed to be instantaneous events. This model leads already 
to possible super diffusion of time-fractional type. The authors are, however, 
looking for a super diffusive regime of Levy type (space-fractional diffusion). 
They then modify their model. The jump PDF and the waiting time PDF of 
the CTRW are now linked by assuming that the particles move during their 
"jump" with final velocity ±V: thus longer jumps are "penalized" by longer 
durations ("Levy walk" rather than "Levy flight", [5]). 

It appears to us that the limitation to two states makes the model un- 
adapted to describing the motion of a physical particle; it is, however, nec- 
essary for the definition of the Levy walk in their model. In the next section 
we show that the assumption about the two states can be dropped, and the 
velocity can be considered as a general random function of time: v(t) [Sees. 
4-5]. 

The model can be further extended by considering a velocity field: v(x, t) 
(Sees. 6 and following). But this case leads to additional difficulties, because 
it requires the determination of the Lagrangian velocity correlation, a prob- 
lem that is well known as a basic difficulty of turbulence theory. Here we 
shall consider only the crudest approximation for the latter quantity, but the 
application of more precise methods that were recently developed in this field 
can be envisaged. 

2 CTRW and Fractional Calculus 

In the forthcoming text, the abbreviation "PDF" denotes a "probability 
distribution function"; the symbol may equivalently be translated as "prob- 
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ability density function" . 

In the Continuous Time Random Walk model one considers a particle 
which at time zero makes an instantaneous jump from position x to x + 
x, then waits at the new position during time t, then jumps again to a 
new position, where it waits, etc. Both x and t are random quantities. 
In all forthcoming calculations, we assume that the random processes are 
symmetrical (centred), i.e., (x(t)) = 0, Vt; hence [see Eq. (1)], 5x(£) = x(t). 
In its standard form, a CTRW is defined by two functions: 

/(x): the PDF of a jump defined by the vector x; Fourier transform: /(k), 
ip(t): the waiting time PDF; Laplace transform: ip(s). 

We only consider in the present work "classical" CTRW's, for which /(x) 
and ip(t) are independent quantities. This means that /(x) is independent 
of t, and ip(t) is independent of x. Processes such as Levy walks [5] or the 
vMSC model [15] are not discussed here. 

A derived function that plays an important role is the memory kernel 
<f>(s) (in Laplace representation): 



4,(8) = 



s ip(s) 
!-4>{s) 



(7) 



This standard CTRW has been solved exactly by [3] in 1965. Let n(x, t) 
be the PDF of finding the walker at position x at time t, given that it 
started certainly at the origin x = at time t — [i.e., n(x, 0) = <5(x) 
and n(x, t) — > as |x| — > oo]. Clearly, from this definition, n(x, t) is the 
propagator (or Green's function) of the CTRW process. It is given, in Fourier- 
Laplace representation, by 



n(k, s) = 



i 



* i-iK*)/(k) 

This is the celebrated Montroll-Weiss equation. An equivalent form is: 

1 



(8) 



n(k, s) 



s + <f>(s) 



/(k) 



(9) 



The latter equation is easily transformed into: 



s n(k, s) - 1 = -<f>(s) l-/(k) n(k,s) 



(10) 
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This is the Fourier-Laplace transform of the integro-differential equation 
obeyed by the propagator: 



d t n(x, t) 



dt! 4>{t - 1') 



-n 



*0 + J«fr-W) 



:n) 



This Montroll-Shlesinger (MS) Master equation is non-local, both in 
space and in time. It shows that the rate of change at point x and time t is 
determined both by the past history, whose influence is measured by 4>(t), 
and by the spatial environment characterized by /(x). 

Although Eqs. (8) and (9) are explicit solutions for the propagator (or 
the resolvent, as the Fourier-Laplace transform of the propagator is usually 
called), the expression of the inverse Fourier-Laplace transform in terms of 
known functions of x and t is, in general, difficult or impossible. 

It has been shown in a series of important recent works that this problem 
can be solved under rather general conditions in the so-called fluid limit, 
which is actually an asymptotic limit of large distances and times [[5], [6], 
[13], [16], [2]]. (It should be noted that, in many interesting cases, the notion 
of "long" is ambiguous when there is no intrinsic time scale!) The fluid 
limit is defined as follows (using dimensionless variables k, s, x, t) [d is the 
dimensionality of the system]: 



f(k) = l-|k| Q ..., 0<a<2, k^O (12) 
^( s ) = 1-/..., < /? < 1, (13) 

In the same limit, the memory function is: 

0(s) = s 1 '? (14) 

This is the same as Eq. (17.12) of [2], describing subdiffusion {(3 < 1) 
or normal diffusion (j3 = 1) [see Sec. 3]. Later in the present work we will 
need an extension to the range 1 < f3 < 2. This extension is not trivial and 
requires a closer discussion, which is developed in Sec. 4. Thus, in the fluid 
limit, the CTRW is entirely defined by the two exponents a and 13. 

The (F-L) equation of evolution takes the form: 

s n(k, s) - 1 = -s 1 -? \k\ a n(k, s). (15) 
Its solution, which is the resolvent of the process, is: 
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s 0-l 

n(k,s) = -5 — — - (16) 
v 7 s^ 3 + |k| Q v 7 

Eq. (15) is, in F-L form, a fractional differential equation (FDE), 

corresponding to: 



d a 

— |k| a — > — — — = D M : Riesz fractional derivative 



9|x| 



s^g(s) - s^git = 0) -> ^£)f #(t) : Caputo fractional derivative (17) 

A primer on fractional calculus, including the definition of these objects 
is given in the Appendix A. Thus, Eq. (15) can be rewritten as: 

n(k,s) -s?- 1 = -|k| a n(k,s), (18) 
which is the F-L transform of the FDE: 

°D?n(x,t) = D? x] n(x,t), (19) 
with the initial condition: 

n(x, t = 0) = <5(x) (20) 

Thus, n(x, t) is the propagator (or Green's function, or fundamental 
solution) of the FDE (19). 

Although the Fourier-Laplace inversion of Eq. (16) is difficult, an impor- 
tant property, i.e., the scaling of the solution in x — £ space can be derived 
without performing explicitly the latter operation. The following argument 
is due to Mainardi [9], and is much more compact than the one presented 
in [17], [18] and [2]. It is based on the well-known properties of the Fourier- 
and Laplace transformations: 

w-lHd- (21 > 

where the double arrow connects a pair of Fourier (Laplace) transforms. 
Applying this property to Eq. (16), we find, after a short calculation: 

1 - /k s\ 1 - fbP /a \ 
n^(ox, bt) — - n a ,p [~,- b )=- n a ,, (— k, sj (22) 

Applying the relations (21) to the last form we find: 
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Thus, finally: 

n ai/3 (ax, bt) = ^ n a ,p (^x, tj . (23) 
Let us assume now that the density profile has the scaling form: 

n atfi (y:,t)=t^ K^J^) (24) 



Under what condition is this form compatible with Eq. (23)? The left 
hand side of the latter yields: 

, . 1 / ax \ 
n a a (ax, bt) = K a g ; 

and the right hand side yields: 

1 / ° x 



,- 7 „ f P7^\ _ 1 ~ / ax \ 



These two expressions are equal, provided that 7 = (3 /a. We thus con- 
clude: 



n 



«A*,t) = t~ Pla Kef* (^) . (25) 

This is the important scaling relation obeyed by the density profile. 
It tells us, in particular, that it depends on space only through the combi- 
nation p = (x/t' 3 /"), called the similarity variable. From this result we 
immediately deduce the scaling of the MSD: 

With the substitution p = x/t' 3 /" we find: 

{x\t))^ = M a ,g t 2 ^, (26) 
identical to Eq. (1), with the following definition of the constant: 

M a , = J dpp 2 K atf) (p). (27) 
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Eq. (26) provides us with a criterion of " strangeness" expressed in terms 
of the two parameters a, (3: 



Diffusion exponent 


Regime 


2(3 < a 


SUBDIFFUSIVE 


2(3 = a 


(NORMAL) DIFFUSIVE 


2(3 > a 


SUPERDIFFUSIVE 



(28) 



The dimensional argument leading to Eq. (26) does not tell us anything 
about the constant M a ^: we do not even know wheter it is finite or infinite! 
The determination of M a ^ requires the explicit form of the solution through 
the function K a ^(p). 

We proceed here in a more direct way [19]. The inverse Laplace transform 
of Eq. (16) is known [7], [9]: 



n^fot) = Ep(-\k\ a 1?), 



(29) 



where Ep{z) is the Mittag-Leffier function [20], a natural (fractional!) 
generalization of the exponential: 



E„(z) 



n=0 



T((3n+1)' 



(30) 



The first terms of the expansion of the Fourier- density profile are thus: 



n a)/? (k,t) = 1 + 



-\k\ a t? 



\k\ 2a t 2lS 



+ 



(31) 



r(a + l) r(2a + l) 

The MSD can be determined directly from the Fourier- density profile by 
the well-known formula: 

<^ = -{s-s^«)L- (32) 

When applied to the Mittag-Leffler function, this yields (in d dimensions) : 



a(a-2 + d)^- — — - 2a{2a -2 + d y-j^ . - 4 



r(/3+i) 



r(2/3 + i) 



Two cases are to be distinguished: 
a) a — 2. 



k=0 

(33) 
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In this case, the first term has a finite positive value in k = 0; all the 
other terms in the series vanish in this point. Hence: 



W))^ = fWTT) a = 2 > W(3 - (34) 

This completes Eq. (26) [Note: {2(3 /a) = (3}. 
b) < a < 2. 

The first term in the right hand side of Eq. (34) contains the factor |/c| a_2 
which diverges in k = 0. The remaining terms are either or oo in this 
range of a. Thus, the constant in Eq. (26) is: 

M a>f3 = oo, < a < 2, V/3. (35) 

This means that the density profiles of the CTRW in the present range 
of a have long tails, whatever the value of (3. Restricting the concept of su- 
perdiffusion to the sole discussion of Eq. (26) would put all the processes 
with < a < 2 on the same rank. The distinction among them requires 
a deeper discussion of the density profile. These questions will be further 
discussed in the next Section. 



3 Classification of Strange Transport Regimes 

The very clear papers of Mainardi and coworkers [8], [9], [10] contain an 
exhaustive study of the Green's function of the FDE (19) in one dimension, 
d — 1. In the forthcoming work, we take full advantage of the fact that 
the general solution n Qii g(x, t) (for all a, (3) of Eq. (19) has been obtained 
by these authors.. This solution was also obtained by [5] in a different form 
(in terms of Fox functions). We shall first discuss some particular cases (for 
d = l). 

a = 2, (3 = 1 

This is the classical case of normal diffusion. Eq. (15) becomes in this 
case: 

s n(k, s) — 1 = — \k\ 2 n(k,s), (36) 

which is simply the F-L transform of the ordinary diffusion equation (with 
a diffusion coefficient D — 1): 

d t n(x,t) = V 2 n(x,t) (37) 
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It is well known that in this case the MSD is: 

(x 2 (t)) = 2t, (38) 
in agreement with Eq. (26). The solution of Eq. (36) is: 

n{hs) = -^- 2 (39) 
The propagator associated with Eq. (37) is the Gaussian wave packet: 

n(x ' <) = 2R^ exp B) (40) 

Note that Eq. (40) is a special case of the scaling relation (25): 

n 2 , 1 (x,t)=rW F{xr l l 2 ) (41) 

The Gaussian scaling function F{p) in the present case is a consequence 
of the Central Limit theorem of probability theory. " Most" physical systems 
behave in this normal diffusive way. But it is the strange behavior, which 
deviates from this Gaussian form, that will be the main object of our study. 

a = 2, < (3 < 1 

This case was studied in detail in [2], and more completely in [8] and 
[9]. It was called the Standard Long Tail CTRW (SLT-CTRW) in [2], and 
time-fractional subdiffusion in [9]. The main result in this connection may 
be stated as follows. Eq. (26) shows that the SLT-CTRW is a subdiffusive 
process (because 2(3 < a), for all values of (3 in the range (0, 1). The 
propagator is non-Gaussian, and has a long tail of stretched exponential form: 

n 2 ,,(x, t) = A % , t^-D/CM) exp (-b, p 2 ^) , p » 1 (42) 

The values of the constants A a ^ and bp [9] are irrelevant here. 

Eq. (42) has again a scaling form, with a (rather complicated) finite 
constant A 2t p. This function has a single maximum in x — 0, and has a 
finite MSD. ' 

0< a < 2, (3 = 1 
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This is the case of space-fractional diffusion of Mainardi 
12.3 of [2]]. The corresponding resolvent is: 

n(k, s) = TTi - - 

v ' s+ \k\ a 

Its inverse Laplace transform is: 

n{k,t) =exp{-\k\ a t). (44) 

This is the Fourier transform of a Levy stable distribution ([2], Appendix 
D). It is well known that the inverse Fourier transform of these distributions is 
in general not expressible in terms of simple functions, except for a few values 
of a. The general propagator (for arbitrary a) can, however, be expressed 
as a Mainardi function, defined by a convergent power series [see Eq. (77) 
below]. The most striking property of the Levy distributions is their "fat" 
long tail, i.e., an algebraic power law decay for p >> 1: 

F(p)~-L;, P»l- (45) 

This law of decay is much slower than the stretched exponential decay 
of Eq. (42). As a result, the MSD of all the Levy distributions is infinite, 
as shown in Sec. 2. The three cases discussed above are illustrated in the 
following figure by the corresponding Mainardi functions. In these figures 
t — 1, hence p = x. 

Red (solid) line: subdiffusive SLT: a = 2, /? = 1/2 

Blue (dotted) line: superdiffusive Levy flight: a = 3/2, (3 = 1. 

Black (dashed) line: normal diffusive Gaussian process: a — 2, (5 — 1 

We represent in Fig. 1 the propagators for these three cases. 

The asymptotic behavior is most strikingly seen in log-log representation 
(Fig. 2). The upper reference straight line has a slope [—(1 + a) — —2.5]. 

We thus see that there are two types of long tails: 

• The stretched exponential tail, characteristic, typically, of the subdiffu- 
sive SLT 

• The algebraic tail, characteristic of the superdiffusive Levy process. 

Note that both are "fat" tails, i.e. the decay at large x is slower than the 
corresponding diffusive Gaussian. The second type is the most slowly decay- 
ing one. Warning: it should not be concluded that all subdiffusive processes 
have exponential tails, whereas all superdiffusive ones have algebraic tails!! 
(see below). 



[9] [see also Sec. 

(43) 
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Figure 1: Propagators (Green functions) G(x) = n(x, t = 1) for three cases of 
fractional diffusion. Dotted blue: Levy super diffusion (a = 1.5, (3 = 1); Solid 
red: SLT subdiffusion (a = 2, (3 = 0.5), Dashed black: Gaussian diffusion 
(a = 2, 13 = 1). 
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a = 2, 1 < (3 < 2 



This case corresponds to time-fractional superdiffusion. This very 
interesting case will be treated in detail in Sec. 5. 



0< a < 2, < (3 < 2 



This is the generic, mixed spatial and temporal fractional diffusion. In 
this case there is a competition between super- and subdiffusion, expressed 
by the ratio 2/3 /a [see Eq. (26)]. In the work of del-Castillo Negrete et al 
[16], a model with a = 3/4, (3 = 0.5 has been successfully compared with 
numerical simulations of a set of plasmadynamical equations. In a recent 
paper [21], additional values of a < 2 were explored. 

In the present work we do not consider the case of space-fractional (Levy) 
superdiffusion. We thus fix in all subsequent sections: 

a = 2 (46) 



4 Random Time-dependent Velocity 

We consider a particle moving in one dimension; its position X(T) at time T 
is determined " semi-dynamically" by a V-Langevin equation ([2], Ch. 11): 2 

^ = V(T), (47) 

where the random velocity V(T) (depending solely on time) is defined 
statistically by specifying the first two moments of the ensemble of its real- 
izations, supposed to be Gaussian, homogeneous and stationary: 

(V(T)) = 0, (V(0) V(T)) = V 2 T(T). (48) 

The function T(T) need not yet be specified explicitly, except by requiring 
T(0) = 1. Vo, the initial velocity rms is supposed to be a given parameter. 

Let f(X,T) be the (fluctuating) distribution function of the particle. It 
is represented as usual as: 



2 Capital notations X,T,V, k, S denote dimensional quantities. An exception is made 
for the wave vector k, because the symbol K will be used later for the Kubo number, Eq. 
(91). 
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f(X,T) = n(X,T) + 6f(X,T) (49) 

where the ensemble-average n(X,T) = (f(X,T)) is called the density 
profile, and Sf(X,T) is the fluctuation. The distribution function obeys the 
Hybrid Kinetic Equation (HKE): 



d T f(X,T) + V(T) d x f(X,T) = 0. 



(50) 



The characteristic equation of this first-order partial differential equation 
coincides with the V-Langevin equation (47). This equation is split, as usual 
[2], into two components by the averaging operation: 



d T n{X,T) +d x (V(T) 5f(X,T)) = 



(51) 



drSf(X,T)+V(T) d x 5f(X,T)+V(T) d x n(X,T)-d x (V(T) 5f(X,T) = 

(52) 

The latter equation is solved by using the propagator: 



g(X,T\X',T') = 5 



X-X'-f dT x V(Tx) 

JT' 



(53) 



Substituting its solution into (51), we find: 



d T n(X,T) = d\ £dT x (v^ViTjn X - £ dT 2 V (T 2 ) , T\ ^ 



+d x ( V(T) 5f 



X 



Jo 







(54) 



This is the exact solution of the problem. It is expressed by an integro- 
differential equation that is non-local both in time (non-Markovian) and in 
space (non-local). We do not repeat here the arguments justifying the lo- 
cal approximation, by which the fluctuating displacements f£ dT 2 V(T 2 ), 

J Q T dTi V(Ti) in the arguments of n and of 5f are neglected, and by which 
the second term of the right hand side, containing the initial fluctuation, is 
also neglected [2]. The result is an (approximate) local, but non-markovian 
diffusion equation: 



d T n(X,T)=d 2 x f dT, {VWViTj) n(X,7\). 
Jo 



(55) 



15 



Using Eq. (48), this becomes: 

d T n(X, T) = d 2 x f T dT x V 2 T(T - 7\) n(X, 7\). (56) 
Jo 

Consider, on the other hand, a particle performing a one-dimensional 
CTRW characterized, in the fluid limit, by a Gaussian jump PDF: 

/(«) = I -a 2 k 2 + ... (57) 

where a is a characteristic length and k is the dimensional wave vector. 
The form of the waiting time PDF is not specified at this stage. 

The density profile 3 (in F-L representation) obeys Eq. (10) which, written 
with a dimensional Laplace variable S, combined with (57) reduces to: 

S fi( K , S) - 1 = -f(S) a 2 K 2 fi( K , S) (58) 
Its inverse F-L transform is: 

<9 T n(X, T) = a 2 d 2 x [ dT x <f>(T - T\) n(X, T ± ) (59) 
Jo 

This equation constitutes the main result of the first part of our work. 
Comparing it with Eq. (56), we see that the density profile produced (in 
the local approximation) by the semi-dynamical V-Langevin equation with 
a time- dependent V(t), and the density profile produced by a CTRW with 
a — 2, - or, more specifically, by the fractional diffusion equation (58) - 
obey formally identical equations of evolution. This statement will be made 
more precise in the discussion of Sec. 5. The quantitative correspondence is 
obtained by identifying the velocity correlation function of the V-Langevin 
process with the memory function: 



V 2 T{T)^a 2 0(T), (60) 
i.e., in F-L representation 

f{S) 4>{S). (61) 



3 The density profile considered in all forthcoming developments is the solution of the 
evolution equation with the initial value: n(X, 0) = S(X). In other words, n(X, T) is the 
propagator of the equation. We shall continue, however, to call this function the " density 
profile" , which has a more physical connotation. 
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We note that Eqs. (56) and (58) were also obtained in Ref. [14] in a quite 
different manner, which is only valid for the 2-state CTRW. The relation (61) 
is much more general and was obtained in a simpler direct way. 

5 Algebraic Velocity Correlation 

We now make the problem explicit, by adopting the same form of the velocity 
correlation as in Ref. [14]: 

V 2 T(T) = A o < 7 < 1, (62) 

where A is a quantity having dimension [length] 2 [time]~ 2+7 . We then 
associate with the V-Langevin equation (for d — 1) (47) [or the Hybrid 
kinetic equation (50)] a CTRW with memory function given, in the fluid 
limit, by Eq. (60), leading to the fractional diffusion equation (FDE) (58): 

S S) - 1 = -A k 2 S 7 " 1 fi( K , S), . (63) 
This relation agrees with Eq. (33) of [14]. It is trivially equivalent to: 

S 2 ^ E(k, S) - S 1 ^ = -A k 2 S) (64) 
The F-L density profile is obtained from Eq. (63) in the form: 

= sJtIak* ' 0<7<l < 65 ) 

Comparing this result with Eq. (16), we immediately see that the two 
expressions become identical by taking the following value for the waiting 
time exponent (3: 

P = 2- r (66) 

Thus, the process corresponds to time-fractional superdiffusion as 

defined in Sec. 3: 

1< (3 < 2. (67) 

At this point there appears a subtle difficulty (the following argument 
results from a remark of J. Klafter [22]). The FDE (63) is formally derived 
as the fluid limit of a CTRW with $(S) = 1 - S?, and f3 in the range (66): 
this is an extension of the range defined in Eq. (13). 
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The knowledge of the Laplace transform ip(S) allows us to calculate the 
average waiting time between two jumps in the CTRW, by using a well-known 
formula: 

(T) = - d ^ {S) 



dS 



= PS< I - 1 \ 8=0 (68) 

s=o 

For < (3 < 1 the average waiting time of the CTRW is infinite. This 
denotes a waiting time PDF with a long tail, leading to strange transport, 
i.e., subdiffusion (Sec. 3 and [2]). For (5 — 1 the average waiting time is 
finite, corresponding to normal diffusion. 

For 1 < (5 < 2, the average waiting time is zero, thus leading to a para- 
doxical situation. Recalling that T is a semi-definite positive variable, it 
follows that ip{T) ~ S(T). Thus, any random walker would remain trapped 
forever at its starting point, in clear contradiction with superdiffusion, as 
predicted by the diffusion exponent fi = 2j3/a = [3 > I (Sec. 3). The only 
other possibility of obtaining a zero average waiting time would be to have 
ip(T) negative in part of the range (0, oo), but then ip(T) does not qualify 
for a probability distribution function. 

The conclusion of this argument is the following. The function u(k, S) 
defined by Eq. ( 65) is, indeed, the resolvent of the fractional differential 
equation (FDE) (63), but the latter cannot be derived as the fluid limit of 
a CTRW. 

We are now faced with a new subtlety 4 . The function u(k, S) is the 
solution of the linear algebraic equation (64). But, is the latter equation the 
F-L transform of a fractional diffusion equation? One is tempted to associate 
Eq. (64) with the FDE in X — T representation: 

) 2 

Vi 



°D 2 T ^n{X, T)=A Df x] n(X, T), < 7 < 1, (69) 



with the initial condition: 



n(X,0) = 6(X). (70) 

This fractional equation is an "intermediate" between a diffusion equation 
(7 = 1) and a wave equation (7 = 0) [23], [8], [10]. But, using the general 
formula for the Laplace transform of the Caputo fractional derivative [7], [10] 
and the Appendix of the present paper, the F-L transform of Eq. (69) is: 

4 We are indebted to F. Mainardi for calling our attention on this point and on Ref. 
[10]. 
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S 2 ~~< n(«, S)-S 1 -^n{K, T = 0)-S^n T (K, T = 0) = -Ak 2 u(k, S), < 7 < 1 

(71) 

It thus appears that the general solution of the Cauchy problem for Eq. 
(69) requires two initial conditions (like the wave equation): n(X, 0) and 
tit(X, 0) = [dn(X,T)/dT]T=o- As a result, the general solution is expressed 
in terms of two propagators; indeed, the solution of Eq. (71) is: 

. We refer to Ref. [10] for a complete solution and discussion of this 
problem. Here we shall only consider the simpler special case (also considered 
in [9]): 

n(X,T= 0) =8(X), n T (X,T = 0) = (73) 



n(«, T = 0) = 1, n T (K, T = 0) = (74) 

In this case, Eq. (71) reduces to Eq. (64), and Eq. (65) represents the 
F-L propagator of Eq. (69) equipped with the initial conditions (73). The 
explicit solution for the propagator of (69) (i.e., the F-L inversion of Eq. 
(65) was obtained analytically by Mainardi [23] , [9] (and also by Metzler and 
Klafter [5] in terms of Fox H-functions): 



n(X,T;i) = — =^ f — X ), 0<7<1 (75) 

This has precisely the scaling form (25), with the similarity variable: 

p= VAP-'nm (76) 

The special function Mi_( 7 / 2 )(p) is called the Mainardi function: it is 
defined for any order v G (0, 1) as: 



/ \ k 

^M=g t , r[ jr + (!-,)) ■ °<"< i < 77 > 

This is a special case of a Wright function, defined as follows ([20]): 



«°-ftri = g «T^+fl - a " s>0 ' (78) 



k=0 
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The series (77) is convergent, but the convergence is very slow. In prac- 
tice, for the numerical calculations displayed below, we take 100 terms in the 
series. Mainardi et al 2001 [9] also showed that the asymptotic behavior 
of the propagator is given by a stretched exponential, which for the value 
(3 = 2 - 7 is: 

n(X,T ;i ) = 2VX *i-(y/2) B pa exp( "^ C) ' P » (79) 

with: 



B 



a = 



b = 



[2tt 7 2 (2 ~^ (2- 7 ) 2 ^- 1 )/t]- 1/2 , 
1-7 



7 



7 
2 

7 



(2-7) 



2-7H 1/7 



< 7 < 1 



(80) 



The important point here is that the present density profile decays asymp- 
totically according to a stretched exponential law, with exponent 2/7 > 
2. This implies that the propagator decays asymptotically faster than the 
Gaussian (contrary to the subdiffusive case < (3 < 1, discussed in Sec. 3. 
Mainardi calls this feature a " thin tail" . It thus appears clearly that superdif- 
fusion is not necessarily associated with a long ("fat") tail of the propagator. 

Another important property of the propagator (75) [i.e., of the Mainardi 
function M v {p) for 1/2 < v < 1, i.e., 1 < (3 < 2] is its possessing for 
T > two symmetrical maxima that move away from the origin, while 
becoming wider. As pointed out above, the fractional equation (69), written 
in a simpler notation {d^/dT 13 = qDj) as: 



n(X, T) = AV 2 n{X, T), l<{3<2 



(81) 



is a kind of interpolation between a diffusion equation {(3 — 1) and a 
wave equation {(3 — 2). Thus, the solution represents a wave propagating in 
both directions away from the origin, and damping out on the way 5 . We 
illustrate these properties by the plot of Fig. 3, corresponding to 7 = 0.5, 
hence (3 = 1.5. For T — > 0, the propagator tends toward n{X, 0) = S{X), as 
it should (70). 



5 In two dimensions, this is a very familiar picture: it represents a wave propagating 
radially on the surface of a lake, when a pebble is thrown into the water. 
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Figure 3: Propagator n(X,T) at four different times. 7 = 0.5, A = l. 

We show in Fig. 4 the tail of the propagator (for T = 1) compared to 
the tail of the corresponding Gaussian. This situation is to be contrasted 
with Figs. 1 and 2: both the time-fractional subdiffusive propagator and the 
space- fractional superdiffusive (Levy) propagator have "fat" tails. 

Figs. 3 and 4 allow us to understand how superdiffusion is possible in 
spite of the thin exponential tail. The dispersion of the matter is produced 
not only by the broadening of the initial distribution (as in ordinary diffusion) 
but also by the symmetric wave-like outward motion of the maximum. 



The results of the present section are not really new: All these results 
were also obtained in Ref. [14], using, however, the 2-state model. We showed 
here that the results are easily extended to an arbitrary fluctuating velocity 
depending solely on time. Our main result can be summarized as follows. 

The motion of a particle in one dimension, with a random velocity depend- 
ing solely on time, obeying a semi- dynamical V-Langevin equation (V-LE) 
(47), can be equivalently described by a time-fractional differential equation 
(FDE). There is, however, no CTRW whose fluid limit is the FDE (74). 
The V-LE is characterized by a velocity correlation depending (for all posi- 
tive times!) algebraically on time, T(T) ~ T~ 7 , with an exponent 7 G (0, 1). 
The equivalent FDE is characterized by a spatial exponent a = 2 and a 
waiting time exponent (3 = 2 — 7; thus 1 < (3 < 2. This FDE leads to 
time-fractional superdiffusion, with diffusion exponent \i — 2 — 7 > 1. The 
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Figure 4: Tail of the propagator for 7 = 0.5, compared to the corresponding 
Gaussian packet. T — 1. 

density profile has two maxima that are moving symmetrically away from 
each other, while broadening with increasing time. The asymptotic decay is 
exponential, decreasing faster than a Gaussian in space ("thin tail"). 

The model V-LE considered here, though allowing an exact analytical so- 
lution, has a physical drawback. The algebraic law of the velocity correlation 
T(T) = T~ 7 cannot be valid for short times. For an obvious physical reason, 
lim T{T) = 1 [as required by Eq. (48)], instead of diverging as in the present 

model. The model should certainly be improved by using a correlation func- 
tion that satisfies this constraint. 

We now generalize our model by considering, instead of a velocity depend- 
ing on a single time variable, a velocity field in two dimensions, for which the 
velocity depends on both time and space V(X, T). 

6 Random Velocity Field with Algebraic Time- 
correlation 

We now return to Eq. (3), and consider a true velocity field, i.e., a function 
both of space and time. We immediately note that this generalization intro- 
duces nothing new in one dimension. We recall the important condition of 
zero divergence (3) which ensures the correct physical status of the hybrid 
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kinetic equation, i.e., the equivalence of Eqs. (4) and (5). Clearly, the con- 
straint (3) reduces, in one dimension, to dxV(X, T) = 0, hence the velocity 
can only depend on time, and we are brought back to the situation of Sec. 4. 
In order to study a non-trivial problem, we must consider d>2. Specifically, 
we shall consider in the rest of the present work two-dimensional systems, 
d = 2. As will be shown below, such systems are actually of much greater 
physical interest than one-dimensional ones. 

We thus consider a particle whose position is specified by the vector X = 
(X,Y) and the velocity V = (V x ,V y ). The corresponding (vectorial) V- 
Langevin equation is: 

^P = V[X(T),T], X(0) = 0. (82) 

As a result of the previous discussion, we may associate with this V- 
Langevin equation a hybrid kinetic equation, provided that the velocity field 
has zero divergence: 

V • V(X,T) = d x V x (X,T)+dYV y (X,T) = 0. (83) 

Contrary to the one-dimensional case, this condition can be realized with 
a non-trivial velocity field. The associated hybrid kinetic equation possesses 
both required properties: it conserves the normalization of / and its charac- 
teristic equations are Eqs. (82). Thus, the two following forms are equivalent: 

d T f(X, T) + V- [V(X, T) /(X, T)\ = 0, (84) 

d T f(X, T) + V(X, T) • V/(X, T) = (85) 

We note that the condition of zero divergence of a two-dimensional veloc- 
ity field is realized, in particular, by the intensely studied model of a charged 
particle moving perpendicularly to a strong constant magnetic field B, in 
presence of a fluctuating electrostatic field (see, e.g. [2], where many other 
references are given]. The coordinates of its guiding centre are animated by 
the well-known electrostatic drift velocity (in Gaussian units): 



dX(T) c d®[X,Y,T\ 



dT B dY 

dY(T) c d$[X,Y,T\ 
dT B dX 



X=X(T) 

(86) 

X=X(T) 
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where $(X,T) is the (fluctuating) scalar potential. The statistical prop- 
erties of the velocity field are thus derived from the statistical definition of 
the potential. We assume that the latter is a Gaussian centred, homoge- 
neous, isotropic and stationary random process, defined by a vanishing first 
moment and by the following factorized form of the second moment: 



($(X,T)> = 0, 
($(0,0) $(X,T)) = £(X) T(T). (87) 

The second relation defines the Eulerian potential correlation, i.e., 
the correlation between the values of the potential at the origin and at a fixed 
point X at time T. In the present work we consider only the case where the 
spatial part of the correlation has a Gaussian form: 



£(X) = exp 



X 2 + Y' 



where A is the so-called correlation length. The temporal part T(T) will 
not yet be specified explicitly. It should, however, satisfy the condition: 
T(0) = 1, and is supposed to depend only on the absolute value of T. 

It is convenient to introduce from here on dimensionless quantities (de- 
noted by lower case letters), defined as follows: 

*=f, V = j, t=^r <P = e- l §. (89) 

Here £ is a measure of the intensity of the potential fluctuations, e.g., 
the root mean square value of the potential fluctuations. T c is an intrinsic 
characteristic time, to be defined later [Sec. 8]. The V-Langevin equations 
of motion (86) reduce to: 



dx.it) ^ r . . , 
-£ L = K v[x(t),t], 

^(x, t ) = -«, *V(M) = « (90) 
The dimensionless constant K is called the Kubo number: 

(91) 

Given the form of Eqs. (90) it is even more convenient to use a further 
reduction of the time variable: 
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6 = Kt. (92) 
The V-Langevin equations then reduce to: 

^ = v[x(6>),0], x(0) = 0. (93) 

The relevant dimensionless Eulerian correlations are now defined as fol- 
lows: 

(0(0,0) 0(x,0)> = £(x)T(rt), 
(v j (O,O)v n (x,0)) = £j n {x) T(K" l 6), (94) 

with (88): 



d 2 

£xx(x,y) = -Tj-^£(x,y) = (1 -y 2 )£(x,y), 
d 2 

£ yy (x,y) = —j^£(x,y) = (1 -x 2 )£(x,y), 

d 2 

£ xy = £yx = T^Q-£(x,y) = xy£(x,y). (95) 

A first, well-known general result is obtained from the V-Langevin equa- 
tion (93): the mean square deviation (MSD) in the x-direction is expressed 
as follows: 

f 9 

(x 2 (6)=2 / dr{e-r)C xx {r). (96) 
Jo 

Here C xx (6) is the Lagrangian velocity correlation, calculated along 
the trajectory x(6>) of the particle, i.e., using the solution of Eq. (93): 

C xx (9) = (v x (O,0)vM9),9]). (97) 

The MSD is related to the dimensionless running diffusion coefficient 
T> x (6) in the x-direction by the well-known Einstein relation: 6 



6 In the general case, the factor 1/2 in Eq. (98) is to be replaced by l/2d. In the present 
situation, although the system is two-dimensional, we are looking for the diffusion in the 
single direction x, hence we must take d = 1. If we calculated (r 2 {9)) , we should take 
d = 2. 
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2 dO 

We now consider the evolution of the density profile n(x, 9) = (/(x,#). 
The distribution function is governed by the hybrid kinetic equation (84), 
or equivalently, (85). The latter is treated by a straightforward generaliza- 
tion of the reasoning leading from Eq. (50) to (55); it will not be repeated 
here. The final result differs from the latter equation in that the simple one- 
dimensional velocity correlation (v{9\) v(9)) is replaced by the Lagrangian 
velocity correlation tensor (vj[x.(9i), 9i] i>fc[x(0), 9}) = Cjk(9 — 9\): 

<%n(x, 9) = Vj [ dr C jk (r) V k ra(x, 9-r). (99) 
Jo 

Clearly, the relation (96) can be obtained directly from Eq. (99). 

7 The Corrsin Approximation 

The study of transport is now reduced to the (difficult) determination of the 
Lagrangian velocity correlation. More precisely, we want to find a relation 
between the (given) Eulerian velocity correlation and the Lagrangian one. 
In the present work we will use only the simplest approximation method for 
this problem: the Corrsin approximation [24], [2]. 

We rewrite the exact Lagrangian velocity correlation as follows: 

C jn (0) = { Vj (0,0) v n [x(0),0]) 
= J dx (vj(0, 0) u n (x, 9) <S[x - x(0)]) (100) 

In the Corrsin approximation, the average in the integrand is assumed to 
be factorized as follows: 

C jn (9)*jdx ( Vj (0,0)v n (x,e)) (<S[x-x(0)]) (101) 

The first factor in the integrand is recognized as the Eulerian correlation 
tensor, assumed to be of the form (87), thus: 

C jn (9) » J d-K S jn (x) T{9/K) (5[x - x(0)]> . (102) 

The last factor is evaluated as follows, in the second cumulant approxi- 
mation: 
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(<S[x - x(0)]) = y dk e ik ' x 



-ik-x(0) v 



y dk e ik ' x exp [- | A; 2 (:r 2 (#))] 



2tt (:r 2 (#)) 



cxp 



Thus: 



r(0//o 



dx Sj n {x.) exp 



x 2 + y 2 
2 (:r 2 (#)) 



x 2 + y 2 
2~W)) 



(103) 



(104) 



2tt (:r 2 (fl)) 

We now use the explicit forms (95) for the Eulerian correlation. It is 
immediately seen by symmetry that the non-diagonal components of the 
Lagrangian tensor vanish, and the two diagonal ones are equal to each other. 



with: 



(105) 



ne/K) 

2tt (x 2 (6)} 



dx dy (1 — y 2 ) exp 



- 1 



(x 2 (6)} 



x 2 + y 2 



The integration over x and y is elementary. The result is combined with 
Eq. (96) to yield: 



m = 



T(e/K) 



1 + 2 f°dr (9-t) C(t) 



(106) 



We thus ended with an integral equation for the Lagrangian velocity cor- 
relation in the Corrsin approximation [26], [2]. This equation will be solved 
in Sees. 9 and 10. 

The diagonal character of the Lagrangian velocity correlation introduces a 
significant simplification in the equation for the density profile, which reduces 
to the following non-Markovian diffusion equation: 



d e n(x, 6)= dr £(r) V 2 n(x, 9-t). 
Jo 



(107) 



The diffusion coefficients in the x- and ^-directions are equal to each 
other. We may therefore drop the subscripts in Eq. (98) and write simply: 
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V(0) = f dr £(r). (108) 
Jo 

Whenever this integral converges, its limit for 9 — > oo defines the ordinary 
diffusion constant: V = P(oo). Whenever £> is finite and positive, the regime 
is normal diffusive. 



8 Algebraic Time Correlation: General Prop- 
erties 

We now make the problem more explicit. We consider a class of V-Langevin 
equations (and associated hybrid kinetic equations) for which the temporal 
part of the Eulerian velocity correlation, T(9/K), is a long-tailed function 
that decays algebraically for long times. It is of the general form (62), but we 
take now some additional precautions to ensure a correct behavior at short 
times. We must, in particular, avoid the divergence at T — > 0. Physically, any 
correct (dimensional) correlation function should have the property: T(T = 
0) = 1. We thus define a characteristic cut-off time T c as the value of T at 
which the right hand side of Eq. (62) equals 1. For all times < T < T c 
we put T(T) = 1. Note that a correlation of the form (62) has no intrinsic 
time scale; but the adoption of a cut-off introduces a characteristic time T c . 
The latter allows us to complete the definitions of the dimensionless time 
t = T/T c (89), of the Kubo number (91) and of the dimensionless variable 
(92). We thus adopt the final form of the dimensionless temporal Eulerian 
velocity correlation as: 

f 1, < 9 < K 
T(K- 1 0) = l IO Q>K , 0< 7 <! (109) 

I 07 ' - 

For convenience, 7 will be called the Eulerian exponent; by definition 
it will be limited to the range (0, 1). Anticipating the results of the next 
Section, it will be shown that the Lagrangian correlation function C(9) has 
the same general algebraic form, with different values of the exponent and of 
the coefficients: 

f 1, < 9 < 9 L , 

m = {± h >o, e>e L . < 110 > 

In order to ensure a correct connection of the two parts, we must take: 
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e L = h l '\ (in) 

A will be called the Lagrangian exponent: its value depends on the 
Eulerian exponent 7 and on K. It will be determined in the next sections. 

Before going into specific calculations, we wish to derive an important 
property relating the nature of the macroscopic transport to the value of the 
Lagrangian exponent. Combining the definition of the dimensionless running 
diffusion coefficient (108) with the form (110) we obtain: 

V{d) -\CdTl + f 0L drh6-\ 9>9 L . (112) 
The integrals are easily evaluated: 



9, 0<9<9 L 

= < )i+I L r «- ( -«), A#L (113) 

The more interesting long-time part is rewritten as the sum of a constant 
and of a time-dependent term; using also Eq. (Ill) we find: 



V{9) = — - h 1 ^ + — 9 l ~\ 9>9 L , A^l. (114) 

We now see the existence of two sharply separated ranges of the La- 
grangian correlation exponent A. 

A) < A < 1. In this case the second term in the right hand side of (114) 
is a monotonically increasing function of time. Hence, for sufficiently long 
times this term strongly dominates the first constant term, and the running 
diffusion coefficient behaves asymptotically as: 

V(9) ~ 1_A , A < 1. (115) 
The regime is thus (time-fractional) superdiffusive. 

B) A = 1. Eq. (113) is no longer valid, but the integral in (112) can be 
evaluated as: 

V{9) = tf/A + h In (°), A = l. (116) 
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The regime is "weakly superdiffusive" . 

C) A > 1. The behavior is radically different. The second term in the 
right hand side of Eq. (114) is now a monotonically decreasing function of 
time. For long times this term decreases to zero. Hence the running diffusion 
coefficient tends toward a positive constant: 

V{6) ~ V = —— h 1/A > 0, A > 1. (117) 
e^oo A — 1 

The regime is now normal diffusive for all values of A > 1. The asymp- 
totic T> is the ordinary diffusion constant. 

Clearly, the value A = 1 is a very special one: as A reaches 1 from below, 
the regime changes abruptly from a superdiffusive to a normal diffusive one, 
and remains so for all values of A > 1. This phenomenon is illustrated in 
Fig. 5. 




Figure 5: Running diffusion coefficient T>(9, A) obtained from the kinetic 
equation, for various values of A < 1 and A > 1. 

The difference between the normal diffusive regime obtained from the 
kinetic equation and the subdiffusive regime predicted by the CTRW for the 
same value of A > 1 appears clearly in Fig. 6. 

In order to obtain a clearer picture of what happens at this point, we study 
the density profile. More specifically, we investigate whether its evolution can 
be described by an equivalent FDE, as was found in Sec. 2. 

As shown in Sec. 4, there exists a formal equivalence between, on one 
hand, the non-Markovian diffusion equation (107) obtained from the hybrid 
kinetic equation and on the other hand, the Montroll-Shlesinger equation 
for a CTRW in the fluid limit. It requires the identification (60) of the 
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D(e,A) 



Lambda = 1.5 




6 40 



Figure 6: Running diffusion coefficient for A = 1.5, as obtained from the 
kinetic equation (solid red) and as predicted from the CTRW (dotted blue). 

Lagrangian velocity correlation with the memory function of the CTRW (in 
dimensionless variables) : 



m - 0(0), 



:ii8i 



or, in Laplace representation: 



(119) 

The Laplace transform of the Lagrangian correlation (110) is, for long 



times: 



C(s) ~ s 



A-l 



(120) 



Comparing with Eq. (14) we find that the temporal exponent of the 
"formally equivalent" FDE is: 



(3 = 2- A. (121) 

Substituting this value into Eq. (16) we find a result similar to (65), with 
A replacing 7. The resolvent given by the Montroll- Weiss equation is: 



n(k, s) 



,1-A 



S 2 - A + \k 2 



:i22i 



The difference with (65) is that we should not limit consideration to 
A < 1. We consider, instead the various cases discussed above. 
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A) < A < 1. The situation is the same as in Sec. 5. There is indeed a 
complete identity between the kinetic result and the FDE with 2 — A, i.e., 
1 < (3 < 2. As shown, however, in Sec. 5, there is no equivalent CTRW in this 
case. The regime is one of time-fractional superdiffusion. The density profile 
is provided explicitly by the Mainardi function (75) 7 . It has two moving 
maxima, as represented in Fig. 3 and decays asymptotically as a stretched 
exponential. The diffusion exponent \x = 2 — A. All this is in agreement with 
the kinetic result. 

B) A = 1. The resolvent (122) of the FDE reduces to: 

This is the Fourier-Laplace transform of a Gaussian packet, characteristic 
of a normal diffusive regime. On the other hand, the result (116) shows that 
the running diffusion coefficient obtained from the kinetic equation grows, 
slowly but indefinitely in time. The FDE is no longer equivalent to kinetic 
theory: A = 1 appears as a bifurcation point. 

C) 1 < A < 2. The propagator (122) corresponds now to a FDE with 
a — 2, < P < 1. It describes a time-fractional subdiffusion, which can 
be derived as the fluid limit of a CTRW. The corresponding propagator 
has a single maximum at the origin, and decays according to a stretched 
exponential [see Eq. (42)]. This regime has been studied in detail in [2]. 
On the other hand, the kinetic theory predicts a normal diffusive regime 
throughout this range of A [Eqs. (114), (117)]. There appears now a divorce 
between CTRW, FDE and kinetic theory. A = 1 has indeed the property of a 
bifurcation point: as the parameter A increases from 0; the propagator first 
follows the "FDE branch". At A = 1 it leaves this branch and follows the 
completely different " diffusive branch" . 

The kinetic equation for the propagator in this regime is obtained from 
Eqs. (99) and (110): 



r 9 

d e n{x, 9) = 9 L V 2 n{x, 9) + h dr r~ A V 2 n(x, 9-r), 9 > 9 L (124) 

For long enough time, the equation can be Markovianized by neglecting 
the retardation in the right hand side. This approximation is justified by our 



7 Because of the presence of the cut-off in the present problem, only the asymptotic 
form (79) is really relevant here. 
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knowledge that a constant diffusion coefficient exists. The final asymptotic 
equation is an ordinary diffusion equation: 

d e n(x, 9)=V V 2 n(x, 9), (125) 

with D given by Eq. (117). Nevertheless, for any finite time this Marko- 
vianization is not reliable, because of the non-exponential Lagrangian corre- 
lation, hence the density profile is in general non-Gaussian. 

D) A > 2. In this range the divorce between kinetic theory and FDE is 
complete. The latter looses its meaning. Indeed, it would correspond to /3 < 0, 
which would describe a waiting time distribution that is an increasing func- 
tion of time. Alternatively, its Laplace transform diverges: ip(s) ~ s - '^ — > oo 
as s — > 0. On the other hand, the Markovianization of the kinetic equation 
(124) is more and more justified as C(9) decreases faster for large 9. 

9 Lagrangian Velocity Correlation: First Ap- 
proximation 

We now determine the Lagrangian velocity correlation as a function of time 
and of the parameters 7 and K of the input Eulerian velocity correlation, 
by solving the integral equation (106). The procedure used is one of suc- 
cessive numerical iterations C^ n \9), starting with the natural Eulerian trial 
function £(°>(0) = T(9/K). We checked the convergence of the procedure. 
As can be seen in Fig. 7, the second iteration provides already an excellent 
approximation. 

The first very important qualitative result is the following. An input 
Eulerian correlation (109) with an algebraic tail produces a Lagrangian cor- 
relation having the same type of algebraic tail. This is clearly seen in a log-log 
representation of the result, as shown in Fig. 8. We thus confirm the antici- 
pated form (111). Our main task will be the determination of the Lagrangian 
exponent A, of the corresponding "intensity" h and of the cut-off time 9l as 
functions of the input Eulerian exponent 7 and of the Kubo number K. 

A first glance at Fig. 7 shows that, even for relatively small values of K, 
the Lagrangian correlation deviates very quickly from the Eulerian one: this 
is clearly seen in Fig. 8. Thus, the nonlinearity, measured by K, is far from 
producing just a small perturbation. 

The first iteration is not quantitatively accurate (Fig. 7), but it has the 
advantage that it can be evaluated analytically. It is therefore interesting 
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Figure 7: Solution of Eq. (106) for the Lagrangian velocity correlation. 
Red solid: starting (Eulerian) correlation; Blue dotted: First iteration. The 
second, third and fourth iterations are superposed on the green curve. The 
parameters are: K = 0.2, 7 = 0.1. 



log L® 




Figure 8: Same data as in Fig. 7, in log-log representation. 
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Figure 9: First iterate of the Lagrangian velocity correlation C^(6), (solid 
red) compared to the Eulerian correlation T{9) (dotted blue) at short times. 
K = 0.5, 7 = 0.1. 



to study explicitly some of its properties which provide a global qualitative 
picture. The explicit form is (106): 



T(0/K) 



1 + 2 j°dr (6-t) T(t/K) 



(126) 



The integral in the denominator, combined with the definition (109), can 
be written in terms of two other integrals as follows: 



where: 



1^6) = [ dr (6-t) = \6 2 
Jo 



(128) 



r in 

W) = / dr(6-r) — 
Jk ^ 

1 K^9 2 ^ — K6 + — *— K 2 . (129) 



(l- 7 )(2- 7 ) 1-7 2-7 
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Combining these partial results, we finally obtain: 

c {1 \e) = (i + e 2 y 2 , 6<k, 



(130) 



0j 



1 + 



(l-7)(2-7) 



K 1 



2 4 - 7 

(l- 7 ) (2-7) . 



2 ' 



9 > K. 



(131) 



The behavior of this function for short times is very simple: starting from 
1, it decreases parabolically: 



CW(9) r^l- 29 2 , 9<^K 



(132) 



As 9 reaches K, we find, of course, that the two branches (130) and 
(131) join each other at the same point. But the derivative at 9 = K is 
discontinuous. Indeed, consider the two integrals (128) and (129); one easily 
finds: 



d 
d9 



h{0) 



= K, 



e=K 



d_ 

d9 



= 0. 



e=K 



Thus, £ (1) (#) exhibits a kink at 9 = K, as can be seen in Fig. 9. 
The asymptotic behavior for large 9 is more interesting. Here we must 
distinguish two ranges of the Eulerian exponent 7. 

a) 7 < 1. In this case the term proportional to 9 2 ~ 1 dominates all others 
in the denominator of Eq. (131) for large 9. We thus find: 



C (1 \9) 
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^275)4-27 



(l-7) 2 (2- 7 ) 2 

We thus find, indeed, for very small 7, an algebraic tail of the form (110): 



£«(0) 



04-7 



9 > K, 7 < 1, 



with: 



;i-7)(2- 7 ) 



K~ 



(133) 



(134) 



Next, we note that the Lagrangian exponent in this approximation is: 
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A( 1 )( 7 )=4-7, 7<1 (135) 

Although this value is not accurate (as can be seen from Fig. 8), it has 
the important property (which will be confirmed) that it only depends on 7, 
not on K (for very small 7). 

b) 7 > 1. In this case the dominant term at large 9 in the denominator 
of Eq. (131) is the term proportional to 9. The asymptotic behavior is thus 
different: 



C (1 \9) 

hence: 



^ (1) (g)~ fc(1 ^ 7) , e>K, 7 »1, (136) 

where k^(K, 7) is found from this equation. The important point is that 
the Lagrangian exponent takes a different dependence on 7 for larger values 
of this exponent: 

A«( 7 )~2 + 7, 7 >1. (137) 

Thus, A (1) (7), which was a decreasing function, starts increasing with 7 
beyond 7 = 1. We may calculate analytically the exact dependence of A^(7) 
by evaluating the slope of the log-log graph of C^(9) (see Fig. 8): 

A«( 7 ) = log[£ (1) (100, K, 7)] -log[£ (1) (1000, X, 7)] (138) 

The result is shown in Fig. 10. Comparing this type of figures for different 
values of K, one finds a very weak dependence of A on K for 7 close to 1. 

The important consequence of this new "bifurcation" is the fact that, for 
all values of 7 we find: 

A (1) ( 7 ) > 3, V7. (139) 

The first iteration thus predicts, in accordance with the discussion at the 
end of Sec. 8, that the first iterate of the Lagrangian correlation leads to 
a diffusive regime for all values of 7. More than that: the Lagrangian 
exponent is in the range D), where there is no equivalence between kinetic 
theory and CTRW or FDE. This statement, however, has to be confirmed 
by the final solution. 
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Figure 10: Lagrangian exponent in the first iteration, as a function of the Eu- 
lerian exponent 7. K — 0.5. The left (right) dotted straight lines correspond 
to Eqs. ( 135) [resp. (137)]. 

10 Lagrangian Velocity Correlation. Final Nu- 
merical Calculation. 

As shown in the previous section, the second iterate of Eq. (106) yields 
already an excellent approximation, while not requiring long numerical cal- 
culations (we therefore omit the superscript ( 2 )).The short-time (9 < K) part 
can even be calculated analytically. We thus have: 

C{9) =- — — , 9 < K, (140) 

v ' (1 + 6»arctan6>) 2 v ; 

C(9) = — — ^ 2, 9 > K. (141) 



0' 



1 + K arctan K + 2j R dr (9 - r) £« (r) 



Eq. (131) is substituted in the integral, and the expression is evaluated 
numerically. 

For short times, the Lagrangian correlation behaves as: 

C{9) ~ (1 + 9 2 )- 2 ~ 1 - 29 2 , 9 <^K, (142) 

i.e., exactly as in the first approximation (132). As can be seen in Fig. 
11, the kink which appeared at 9 = K in the first approximation is still 
present in the final expression. 
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Figure 11: Second iterate of the Lagrangian velocity correlation & 2 >(9) = 
C(9), (solid red) compared to the first iterate C^ l \6) (dashed green) and to 
the Eulerian correlation T{6) (dotted blue) at short times. K = 0.5, 7 = 0.1. 

The long-time behavior is, of course, more interesting. The first important 
point, that is seen in Fig. 8, is that the Lagrangian velocity correlation 
behaves asymptotically as a power law, (110), thus confirming our previous 
anticipation. We now study its parameters. 

The Lagrangian exponent A is in all cases smaller than the one given by 
the first iterate, but much larger than the Eulerian exponent 7, as seen in 
Fig. 8. In Fig. 12 we plotted log C(6) vs log (9 for given 7 = 0.1 and three 
values of K — 0.1, 0.5,0.7. For all these relatively small values of K 8 the 
exponent A is strictly the same: A(7 = 0.1) = 2.10, \/K. 

Next, we measured the dependence of A on the Eulerian exponent. The 
result is a very simple linear dependence, shown in Fig. 13. This dependence 
is very accurately fitted by the formula: 

A( 7 ) = 2 + 7, 0<7<1 (143) 

We note already here the following very important point: the Lagrangian 
exponent A (7) > 2, V7. Hence the process described by the V-Langevin 
equation (93), with the Eulerian velocity correlation (110) in the Corrsin 
approximation, leads to a diffusive regime for all values of 7 > and of 
K. This point will be further discussed below. 

8 We know that the Corrsin approximation is not valid for large Kubo numbers. 
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Figure 12: Lagrangian velocity correlation vs. 6, in log-log representation. 
The curves correspond to K = 0.1 (solid red), K = 0.5 (dotted blue) and 
K = 0.7 (dashed green). 7 = 0.1. 




Figure 13: Lagrangian exponent A vs. Eulerian exponent 7. The curves for 
K = 0.1, 0.5, 0.9 are all superposed. 
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We now consider the numerator h of Eq. (110), determined from the 
logarithmic representation of C(6) (Fig. 12). We first note (Fig. 14) that, 
for given K, it is a very slowly varying function of 7: we therefore consider 
it approximately independent of 7. 




Figure 14: The coefficient h(K, 7) vs. the Eulerian exponent 7, for three 
values of K: K = 0.1 (solid red), K = 0.3 (dotted blue), K = 0.5 (dashed 
green) . 

The dependence on K appears to be quadratic, for relatively small values 
of K. It is quite well fitted by the following empiric formula, valid for 7 < 0.5 
(Fig. 15): 

h(K) ~ 0.37+ 1.85 K 2 . (144) 
Finally, the value of the cut-off time is obtained by setting C{Ql) = 1. 

9 L { 1 ,K) = h{Kfl^ (145) 

To sum up, we found a reasonable analytical representation for the La- 
grangian velocity correlation in the Corrsin approximation. A simple feature 
of this representation is the fact that (see Fig. 13 and Fig. 14) each of the 
two output parameters, A, h, is a function of a single input parameter 7 or 
K: 




#<#l(7,^), 
6>6 L ( l7 K) 



(146) 



combined with Eqs. (144) and (145). This analytic approximation is 
compared to the "exact" (numerical) Lagrangian velocity correlation in Fig. 
16. 



41 



Q2 L 




0.2 0.4 



Figure 15: The coefficient h as a function of the Kubo number K. Red 
dashed: 7 = 0.1, blue dotted: 7 = 0.5. The solid black curve is the empirical 
fit ( 146). 
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Figure 16: "Exact" numerical Lagrangian velocity correlation (solid red) 
compared to the analytical approximation (146) (dashed black). 
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11 Running Diffusion Coefficient 



The running diffusion coefficient for an input algebraic Eulerian correlation 
is obtained from Eq. (108). Using the algebraic form (111) this integral was 
calculated analytically in Eqs. (110) and (114). Using also the forms (143) 
- (145) we obtain the result plotted in Fig. 17 for fixed K = 0.3 and three 
values of 7 = 0.1, 0.3, 0.5. As we know that the regime is always diffusive, the 
running diffusion coefficient tends toward a finite positive diffusion constant 
for sufficiently large times (in practice, 9 > 10). One sees that the asymptotic 
diffusion coefficient is a weakly decreasing function of the Eulerian exponent. 




5 10 15 20 

e 

Figure 17: Running diffusion coefficient T>(9) for three values of 7 = 0.1 
(solid red), 7 = 0.3 (dotted blue), 7 = 0.5 (dashed black). K = 0.3. 

Fig. 18 shows the dependence of the diffusion coefficient on the Kubo 
number for fixed 7 = 0.1. The asymptotic diffusion coefficient is a rather 
strongly increasing function of K. 

In most existing works on the Langevin equation the Eulerian velocity 
time-correlation is assumed to decay exponentially (/3 = 1). This produces 
a Lagrangian correlation that also has an exponential asymptotic behavior, 
and thus automatically ensures the convergence of the integral (108) defining 
the diffusion constant. It is interesting to compare the diffusion coefficients 
obtained in the algebraic case and in the exponential case. For the latter we 
use a Eulerian correlation that is tailored in such a way as to resemble the 
algebraic one: we thus start the exponential decay at the cut-off time of the 
algebraic one (Fig. 19): 
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Figure 18: Running diffusion coefficient T>(6) for three values of K = 0.1 
(solid red), K = 0.3 (dotted blue), K = 0.5 (dashed black). 7 = 0.1. 



% xp (e/K) 



exp 



1, 9 < K 
9-K 



K 



6>K ' 



(147) 



The Lagrangian velocity correlation is calculated in the Corrsin approx- 
imation up to the second iteration, by following exactly the same steps as 
in Sees. 9 and 10. The result is shown in Fig. 20. (Note that £ exp (0) = 1, 
hence it is not necessary to introduce a cut-off time 9l as for the algebraic 
case). As expected, C exp (9) decays much more quickly: it has no long-time 
tail as the algebraic one. 

The running diffusion coefficients corresponding to the same parameters 
are compared in Fig. 21. The algebraic correlation produces a running dif- 
fusion coefficient that saturates at a longer time. As a result, the asymptotic 
diffusion constant T> > T> exp . This results from the fact that the area under 
the Lagrangian correlation is larger in the algebraic case (Fig. 20), due to 
the slow decay of the correlation. 



12 Conclusions 

Our main purpose in the present work was to study the relation between 
a "semi-dynamical" description of a system of particles, based on the V- 
Langevin equation and its associated hybrid kinetic equation, and on the 
other hand on a purely stochastic description in terms of a continuous time 
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Figure 19: Temporal part of the Eulerian velocity correlation. Solid red: 
exponential, Eq. (147); dotted blue: algebraic, Eq. (109). K = 0.5, 7 = 0.3. 
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Figure 20: Lagrangian velocity correlation for the case of an algebraic (solid 
red) [K = 0.5, 7 = 0.1] and of an exponential (dotted blue) Eulerian corre- 
lation with the same K. 
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Figure 21: Running diffusion coefficient for the case of an algebraic (solid red) 
[K = 0.5, 7 = O.f] and of an exponential (dotted blue) Eulerian correlation 
with the same K. 

random walk and its limiting form as a fractional diffusion equation. We 
wanted, in particular, to find out under what conditions on the velocity 
correlation function a V-Langevin equation describes a "strange transport 
regime" , which could be equivalently represented by a CTRW or a FDE. 

It was known from previous work [2] that a subdiffusive regime can appear 
in some limiting cases. This occurs, e.g., for a collisional plasma in a strong 
magnetic field, in the limit of infinite correlation length perpendicular to B 
[25], [26] or for electrostatic drift wave turbulence in the limit of infinite 
Kubo number [2]. In these cases, subdiffusion is produced by sticking of the 
particles to the magnetic field lines, or by their trapping in the turbulent 
field. It might be stressed that in the latter diffusive or a subdiffusive 

(K — > oo) regime was found even when the spatial correlation is long-range 
(Lorentzian), but the time correlation is exponential [27]. A superdiffusive 
regime is, however, much more elusive. To the best of our knowledge, a 
V-Langevin equation producing Levy-type superdiffusion was only found for 
a one-dimensional two-state model (velocity = ±V) [14]. The restriction to 
a velocity constant in absolute value was necessary in order to replace the 
Levy flight CTRW by a Levy walk. 

In the present work we showed that, for one-dimensional systems, the re- 
sults of [14] concerning the time-fractional superdiffusion can be extended to 
the more physical case of a general random time-dependent velocity function 
equipped with an algebraic time-correlation function (but with a Gaussian 
space-correlation). We confirm that in this case we find a process that is 
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an "intermediate" between diffusion and wave propagation. This type of su- 
perdiffusive processes, distinct from the Levy-type, has been studied in detail 
in [9] and [10]. 

We then generalized our investigation to the study of a random velocity 
field in two dimensions. In this case one is faced with the difficult problem of 
determining the Lagrangian velocity correlation, a problem that was treated 
within the Corrsin approximation. It was first shown that if the Eulerian 
correlation decays algebraically as Q^ 1 for long times (and its spatial part is 
Gaussian), the Lagrangian correlation has the same property, with a different 
exponent: 9~ A . 

We then established a general criterion for the existence of a superdiffusive 
regime, based on the value of the Lagrangian exponent A. It was shown in 
Sec. 8 that whenever < A < 1, the regime is superdiffusive and can be 
equivalently described by a time- fractional FDE. When A > 1, however, 
the semi-dynamical description predicts a diffusive regime for all A, whereas 
the formally associated FDE predicts subdiffusion, or even looses meaning, 
for A > 2. Thus, the purely stochastic CTRW or FDE descriptions of the 
turbulence are severely limited in this case by the existence of a bifurcation 
at A = 1. 

In the final stage (Sec. 10) the Lagrangian exponent A was determined 
numerically. The result shows that A = 2 + 7 > 2, for all values of 7 and 
of K. Thus, for the present model and with the approximations used in 
its treatment: [local approximation of the non-Markovian diffusion equation 
(54), Corrsin approximation of the Lagrangian velocity correlation (101)]: 

• The regime is always normal diffusive, for all values of the input 
Eulerian parameters 7 and K. There never appears any superdiffusion, 
in spite of the long algebraic tail of the temporal velocity correlation. 
This is to be contrasted with the one-dimensional case of a purely time- 
dependent velocity, where the regime is superdiffusive for all values 
< 7 < 1. 

• The "formally associated" fractional differential equation, has no mean- 
ing. The evolution cannot be described by a purely stochastic process. 

• Although the regime is diffusive, the density profile is not a Gaussian 
packet. Indeed, although the Lagrangian correlation decays sufficiently 
fast in order to ensure the existence of a finite asymptotic diffusion 
constant, the non-Markovian character of the equation (107) cannot be 
neglected for finite time. 

The results obtained here may appear somewhat disappointing, in the 
sense that no time-fractional superdiffusion was found for two-dimensional 
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motion. Nevertheless, we believe that the result of the exhaustive analysis 
performed here is interesting in exhibiting the limits of a purely stochastic 
description. It shows that a long-tailed algebraic temporal correlation is not 
sufficient for producing superdiffusion (for dimensionality d > 1). In order 
to obtain a Levy-type superdiffusive regime, we must presumably act on the 
spatial part of the Eulerian velocity correlation. This problem is, however, 
highly non-trivial. On one hand, in the Corrsin approximation, Eq. (106) is 
invalid, because it is based on the Gaussian form of the spatial Eulerian cor- 
relation. But one should go back even further in the chain of approximations. 
The local approximation leading from Eq. (54) to (55) (in two dimensions) 
is questionable in the presence of long-range spatial correlations. We hope 
to come back to this problem in forthcoming work. 
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Appendix: A Primer on Fractional Calculus 

Fractional Calculus has been known to mathematicians for a very long time. 
The problem appears to have been first raised by L'Hospital and by Leibniz 
in 1695. The first systematic construction appears in the early 19-th century 
and is due to Liouville (1832) (well known to every physicist from the first 
lecture in statistical mechanics), who used this concept in potential theory. 
It was further developed by B. Riemann in 1847. For a long time the concept 
was studied only by mathematicians. It entered physics only in the last half 
of the 20-th century, through the theory of stochastic equations associated 
with the problems of random walks (Sec. 2). 

The search for a natural generalization of the ordinary derivative d n /dx n , 
with n: an integer, is clearly attractive for a mathematician. It is as natural 
as the passage from the integers to rational, and later to real numbers. We 
must find a way to interpolate in a consistent way the concept of differenti- 
ation between two successive integers, n and n + 1. The way of doing this is 
however not obvious. Indeed, several different approaches have been devised 
by mathematicians. We shall only discuss here the one which is most widely 
used in the literature (for a discussion of other approaches, see [7], [28] and 
[11]). One will find detailed expositions of the subject in [8], [13], [7], [29], 
[5], [6], [30]. Excellent tutorials are the papers by [23] and especially [9], [10]. 
A brief, but very clear summary is given in [16]. 

The names of fractional integral, fractional derivative, fractional calculus 
are actually rather improper. The object of these concepts is to extend the 
notions of integration and of differentiation from integer order n to arbitrary 
real order [3, not necessarily a fraction. The name has, however, acquired 
general practice. 

RIEMANN-LIOUVILLE FRACTIONAL INTEGRAL 

Consider first the multiple integral of order n (n: integer) 9 of a real- 
valued function f(x) of the real variable x: 



The two lower subscripts relate to the limits of integration (a is a fixed real 
number); the superscript n denotes the order (multiplicity). This equation 
can be written in the following equivalent form, in terms of a single integral: 

9 In the forthcoming formula we will always use latin letters m,n,r, ... for integer indices, 
and greek letters for real, non-integer indices: a, (3, ... 




(148) 
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••™*>=^f*<^h- (149) 

The equivalence between the two forms is easily demonstrated by partial 

integration. We may put a J° f(x) = f(x) [7]. For n = 1, the formula 
reduces to the single integration: 

a Jlf{x)= f dyf{y) (150) 

J a 

Let us check n = 2: 



«■£/(*) = f X dy(x-y)f(y) 

J a 

/y 
dx 1 f{x 1 ) 

= dx ± dx 2 f(x 2 ). 

J a J a 



-1) dy / dx 2 f(x 2 ) 



In the form (149) the way toward a generalization of the operation of 
integration (149) is obvious: 

aJxf{x) = Jady ^ ^}yy-a ' " > : any positive real number. 

W (151) 
This is the definition of the Riemann - Liouville (RL) fractional 
integral of f(x) of arbitrary real order a. 

Some of the more important properties are given here: 

A) RL fractional integral of a monomial 

° J * ^= rTZ^7T\ x " +a - (152) 
1 {n + a + 1) 

B) Commutative composition rule: 

J£ • „■£ ■ f(x) = oJ v x ■ ■ f(x) = Jf • f(x) (153) 
Proof: 

„j? ■ j: ■ m = ^ f % (x - ^- jT dr « - t)- i+ -/(t) 
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Substitute: £ = r + (x — t)( 



= f^f^ [ dr }(t) £ d( (x - r) [«* - t)]-'+" [(1 + Q(x - r)]-'+' 

- f^m f dT /(r) (I - T) ~ 1+ " + * i 1 iC c_1+ " (1 + 

The ^-integral is well-known is expressed in terms of the well- 
known .B-function: 

jf^.(l- Tr . = fl( ,,,0 = ™. (154) 
hence: 

^■^■ /(l) = i^r dr(«-r)-^/(r) 
= Jff(x) = o^- oJ£ •/(*), 



RIEMANN-LIOUVILLE FRACTIONAL DERIVATIVE 

The fractional derivative is constructed by combining the operations of 
the usual (integer-order) derivative with the fractional integral. 
We define the fractional derivative as follows: 

d rn 

a D a x f{x) = — [ a 4 m - a) f(x)} , a>0, m-Ka<m. (155) 

where a is any positive real number, not an integer, and m is the smallest 
integer larger than a. 

In order to "understand intuitively" this rather unusual definition, we 
may think, very casually, and without any rigour, as follows. We have a solid 
reference point, viz. a rigorously defined fractional object, the RL integral. 
We then start with a RL integral of order m — a, where m is the smallest 
integer which makes m — a > 0, thus ensuring that a Jz™ 1 °^f(x) is a true frac- 
tional integral (151). Next, we "differentiate away" the integer m, by taking 
the m-th ordinary derivative of this object. This looks like "lowering the 
order" of a Jx m ~ a \ thus "transforming" it into a Jx~ a , which could be viewed 
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as a fractional integral of negative order, —a, to be interpreted as a fractional 
derivative of order +a. But, again, this paragraph is just cavalier talk; in 
particular, viewing J~ a as an "inverse" to D a is an ambiguous statement, as 
will be shown below. We now go back to business, and consider Eq. (155) 
as a given definition, and study some of its consequences. 

A basic requirement is that the definition (155) be compatible with the 
ordinary (integer-order) derivative. When a is an integer, a = N, we have 
m — 1 = a = N, thus m — N + 1, and recalling Eq. (150), we have: 



jN+1 jN j px jN 

= ^ M/w] = s I * / w = ^ /<*); 

the operation thus reduces to the ordinary differentiation, as it should. 
Combining Eqs. (151) and (155) we find: 



d m 1 f(y) 

a D" fix) = - — — f x dy r — , a>0, m — 1 < a < m. 

xJK ' dx m V(m-a) Ja y ix-v)^ 1 -™' 

(156) 

This is the definition of the Riemann - Liouville (RL) left-fractional 
derivative of order a. 

We may also define a Riemann-Liouville (RL) right-fractional deriva- 
tive in which it is the upper integration limit b that is fixed: 



d m (-l) m f b 6(v) 

(157) 



We now explore a few properties of the fractional derivatives: we shall 
meet some surprises. These are all due to a fundamental fact. The fractional 
derivative (for any of the forms defined above) of a function f(x) is ex- 
pressed as an integral of f(x) over a given range of x 10 It is therefore 
an intrinsically non-local operator. In other words, the value of the frac- 
tional derivative (for non-integer a) a D%f(x) at the point x depends on the 
values of f taken in the whole range of integration defining this quantity. 
It is quite surprising that, upon varying continuously the order a, this non- 
local character disappears suddenly whenever a reaches any integer value! 

10 For this reason, the fractional derivative is sometimes called " differintegral" . 
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What is not surprising, therefore, is that fractional calculus entered physics 
when researchers began to be interested in evolutionary processes endowed 
with memory ("non-Markovian processes") or in the influence of the spatial 
neighborhood in the phenomena at a given point (x,t). 

a) Linearity: 

D« [A f(x) + g(x)\ = A D° f(x) + » D«g(x) (158) 
The proof follows trivially from the definition of the RL derivative. 

b) RL Fractional derivative of a monomial: 

D° x"= ^ + x»- a , a>0, n>-l. (159) 
1 {/i — a + 1 ) 

Proof. Using the definition (156): 
d m 1 f x 

« D ° x " = ^ rfrzjj I * v <* - vr 1+ - (m 

Let us do this calculation in some detail. Substituting y — > xt we 

find 



/ d m \ 1 r 1 

D« x» = — - x m+ ^- Q — / dt t» (1 - t)- 1+m ' a . 

\dx m J r(m-a)J K ' 

The integral appearing here is expressed in terms of the 5-function, 
(154), hence: 

x \dx m / r(m-a) r(/i + l-a + l)' 

Next, we calculate (setting 7 = /i — a) 

d rn 



x 



m +7 = ( m + 7 )( m _ 1 + 7 )...(1 + 7 ) X 7. 



Noting the identity: 

T(m + 7 + 1) = (m + 7)(m - 1 + 7)...(1 + 7) T((l + 7), 
we combine all these partial results to obtain the final formula (159). 

The limitation to fj, > — 1 in Eq. (159) is due to the factor T(l + ji) 
which has a pole in /i — — 1. Note that for any non-integer value of a the 
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RL derivative qD" x^ 1 has a single zero (for any x 7^ 0) in the allowed region, 
namely for \i — a. — 1 (see below). Indeed, any zero of the right hand side of 
Eq. (159) originates from the factor l/T(p+l — a) which vanishes whenever 
its argument is zero or a negative integer, see Fig. 22; a single one of them 
is larger than —1. 

D« x a ~ l = 0, Va > 0. (161) 




-5 UJ LJ m ^ 1 

"4 "3 -2-10 1 2 3 

Figure 22: Coefficient C(a,/i) = r(/i + l)/T(— a + /i + l) of the RL derivative 
of x M as a function of /i. a — 0.75. 

The formula (159) is deceptively simple, and would correspond with in- 
tuition. Just as the ordinary integral-order derivative {d n / dx n )x m lowers the 
power of x by n, the fractional derivative transforms x M into with a 

prefactor. But a further reflection shows us that this result is not as innocent 
as it appears. 

bl) RL Fractional derivative of a constant: Consider, indeed, the case 
/i = 0, i.e., the 

T(l) 

oD x cx - c ^—-^ x , 

i.e., 

G D« c = ° x~ a ± (162) 
1 1 — a) 
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Thus, the RL fractional derivative of non-integer order a > of a con- 
stant is non-zero! Is this result consistent with the ordinary derivative? The 
answer comes from the properties of the gamma function. Indeed, whenever 
a = n, where n is a non-zero integer, n = 1,2,3, ... the argument of the T 
function is a negative integer, or zero. But we know that 

r(o) = r(-i) = r(-2) = ... = 00 

Hence, for any integer value of a, Eq. (162) reduces, as it should, to: 

££c = 0, n = l,2,... (163) 

Another feature of Eq. (159) seems to contradict our previous statement: 
D" x^ is a perfectly local expression, depending only on the value of x^ at 
x. 11 In order to understand what happens here, we consider the RL derivative 
of x^ for a lower integration limit < a < x instead of zero. Starting the 
calculation as above, we find after the first step: 



D a r M 



dr 



dx r 



x 



T(m - a) J a 



/X 



<% e (i-o 



— l+m— a 



X > a. 



(164) 

There are two important differences with Eq. (160). First, the integral 
depends on x, through its lower limit; hence, the m-th derivative acts on the 
whole product, thus yielding a much more complicated x-dependence. Next, 
the integral with a non-zero lower limit is no longer a 5-function, but rather 
a very complicated hypergeometric function. In order to illustrate our case, 
we consider the simplest non-trivial case that leads to simple calculations: 
/j, — 2, a — |, hence m — 1. Eq. (164) then reduces to: 



D l J 2 x 2 



( iS,, 4l« ,(1 - {r " ! ) 



(165) 



d 

dlc{ x/ r(i/2) Ja/: 

The integral is tabulated; a simple calculation then leads to the result: 



Dl' 2 x 2 



1 



a^x 



X 



\Jx — a 



+ 2x\Jx — a — 



(x - af' 2 



The result reduces to Eq. (159) when a — > 0: 



,Dl' 2 x 2 



r(3) 



3v^ 



x 3/2 = L x 3/2_ 



r(5/2) 



x > a. 



(166) 



(167) 



11 This feature is never mentioned in the literature, as far as we know. 
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The non-local character of the fractional derivative is now luminously 
illustrated. The complicated function appearing in the right hand side of 
Eq. (166) is determined by the choice of the function £ M in Eq. (164), 
hence by its values taken over the whole integration range. In particular, the 
fractional derivative depends explicitly on the lower integration limit a. A 
striking fact is the singularity of the RL fractional derivative when x — > a. 
This singularity is present for any a, arbitrarily close to zero, but disappears 
when a = 0. When x ^> a, a D l J 2 x 2 — > qD]/ 2 x 2 . All this clearly shows that 
the value a = is not a generic one. The situation is illustrated in Fig. 23. 
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Figure 23: Fractional derivative of x 2 : a D 1 J 2 x 2 = D(x, a), for a = 0, 5, 10, 12. 
c) Inversion of the RL fractional derivative. 

Many authors state rather loosely that "the RL fractional derivative is 
the inverse of the RL fractional integral". They even adopt the notation: 
a J" = a Dx a \J\- ft wm be shown now that this statement is, for the least, 
ambiguous. Hence Podlubny's notation should preferably be avoided in order 
to avoid misunderstandings (Mainardi). 

cl) RL Fractional derivative of a fractional integral of same order. 

D: [ J« f(x)} = f(x). (168) 
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Proof: 

For integer a = n > 1: 



= £<^£ *(*-vr'/(») 

For non-integer a, we enclose this number again between its 
nearest neighbor integers, defining the integer A; as: k — 1 < a < k. 

Using now the composition rule (153) for the fractional integrals, 
we note the identity: 

o Jt a ■ [ j: /(*)] = J k x m 

Using now the definition (155) we have: 



d: [ j: /(*)] = — k { Jt a [ j: m}} = ^ J k x nx) = /(*). 

(The result follows from the fact that the last step only involves 
operations of integer order k, i.e., ordinary derivative and integral.) 

From the result (168) one is tempted to conclude that, indeed, "the RL 
fractional derivative is the inverse of the RL fractional integral": the frac- 
tional RL derivative annuls the action of the fractional integral. The situation 
is, however, not so clear-cut. It appears from the following property that the 
double operation D J is not commutativel Indeed: 

c2) Fractional integral of a RL fractional derivative of the same order. 



oJ« [„££ /(*)] = /(*; 



Proof. 



i=i 



x 



a- j 



r(a-j + l) 



k - 1 < a < k. 
(169) 
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oJ: [oD a x f(x)} = -L J X dy (x-y)^ 1 D* f(y) 



i dy — 

T(a) J a 

d 



T x {x ~ v) 



dx \ T(a + 1) Jq 

d r 

dx \ T(a 



oDy f(y) 

dy (x-yr D«f(y) 



A partial integration yields: 



Tx{w)L dy{x ~ y)c 



.! d 

dy 



k-l 



r(« + i) 



x 



d: 



k-l 



dx k ~ 



x=0 



Continuing the partial integrations, and using Eqs. (151) and 
(153), we end up with: 



'"' r fjk-j 



E 

j'=i 



^ r TQ-fe+1 rfe-a 



,=o r(a-j + 2) 



-E 



dx 



oJlf(x) - 

3=1 



d k ~i 
dx k ~ 3 

1 r d k ~ j 



x 



a-j+l 



x=0 T(a-j + 2) 



dx k ~i 



x 



a-j+l 



= m - 



dx k 3 



x 



o T(a-j + 2) 

a-j 



__ o r(a-j + lY 



QED 



It is precisely because of this non-commutativity of integration and dif- 
ferentiation that one cannot state unambiguously that one operation is the 
inverse of the other. 
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d) General composition rules. 

In all subsequent formulae we take the lower limit a = 0. We simplify 
the notations as follows: qJ" = J a , qD™ = D a . The proofs of the following 
results are very similar to the previous ones [7]. 

dl) Composition of two fractional integrals, [see (153)] 

J a J 13 f(x) = J 13 J a f(x) = J a+I3 f(x) (170) 
d2) RL derivative of fractional integral. 

D° J? f(x) = J?- a f(x), 0<a<{3 

D a Jl 3 f(x) = D a ~P f(x), 0<{3<a { > 



d3) Fractional integral of RL derivative. 



a- j 



J a D 13 fix) = D 13 ' fix) — V \Df>-> fix)] n — 

M ' Jy ' p x v JK Jix=0 r(a-j + iy (172) 

< a < (3, n - 1 < (3 < n 



n T a—J 

J a D 13 fix) = J a ~? fix) - V fix)] n — 

M ' Jy ' p x L M Jix=0 T{a-j + l) (173) 

< (3 < a, n - 1 < (3 < n 



d4) Composition of RL derivatives. 



it 

oD« ■ Df> f(x) = D^ f(x) - £ [ oDt j M] x=0 



m 

L>£ ■ D% f(x) = oD^ f(x) - £ [ oDr j f(x 

3=1 

m — 1 < a < m, n 



-a~J + 1) 

x -P-3 



x=0 T(-(3-j + l) 
1 < (3 < n 



(174) 
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Integral transforms of fractional operators. 
A) Fourier transform 

The Fourier transform involves integrals ranging from — oo to +00. 
This implies a number of technical assumptions about the regularity of the 
functions appearing in the theory: we do not discuss them here, assuming 
that all conditions are satisfied for the existence of the objects of the forth- 
coming equations. We begin with the following result: 

e ikx = (ik) a e ikx . (175) 

Proof: 



1 d m r x 

F + (x) = _«,££ e ikx = — / dy(x- y)- l +™-" e «* 

1 (m — a) dx m J _ 00 

1 /jrn poo 

= — e ikx / dt r l+m - a e~ ikt 

r(m-a) dx m J 

, m — 1 < a < m 

where we set x — y = t. We now note that, from the definition 
of m, we have: < m — a < 1. We may thus use the tabulated integral: 

/•OO 

/ dt, c 1 e"^ = (iA;)-^ r( 7 ), 0< 7 <1 
Jo 

Hence: 



1 d m 

= {ik) a e ikx , QED 

This expression is not very convenient in applications. On one hand, the 
appearance of i a is somewhat awkward. But more important, if we wish to 
use this result in relation with the Fourier transform, we need to cover the 
whole range (—00, 00), and not only the left range (—00, x). We therefore 
consider also the right RL fractional derivative (157): 

X D^ e ikx = (-ik) a e ikx . (176) 

Proof: 
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(—'\) m d m f°° 
F_(x) = X D^ e fe = ] — / dy(y-xy 1+m - a e ik * 

1 (m — a) dx m J x 

(—]) m d m r~ x 
Setting £ = — x, we have: 



Thus: 



F_(x) = F;(0 = (-ik) a e^ = (-ik) a e lkx QED 

In order to adapt the previous results to the Fourier transform, in partic- 
ular to cover the whole range (—00, +00), we define a new type of fractional 
derivative, which is totally symmeric. It is justified by the linearity property 
of the RL derivative (158): 

B? x \ /(*) = ^ [ -00^ + X D«J f(x) (17?) 

2 cos „ 

2 

This is called the Riesz fractional derivative. 

The usefulness of this concept appears obviously in the following result, 
expressing the Fourier transform of the Riesz derivative: 

* {D? xl f(x)} = -\k\«f(k) (178) 

where f(k) denotes the Fourier transform of f(x). 
Proof: 

* \P\x\ /(*)} = ^ [ D \*\ I dk'e- ik> * /(F)} 
= -T I J dk' \k'\ a e~ ik ' x f(k') } 



dx e ikx / dk' \k'\ a e~ ikx f(k') 



1 

'27T 

■ |A;| Q /(£;), QSD 
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The result (178) ends our quest of the Fourier transform of the (Riesz) 
fractional derivative of a function, expressed in terms of the Fourier transform 
of that function. 

B) Laplace transform 

The (direct) Laplace transform of a function involves an integral over the 
variable from to oo. It is therefore particularly well adapted to the study 
of functions of time, when the latter is restricted to positive values. The 
two properties of the Laplace transform that are most useful in applications 
(especially for the solution of differential equations are: 

The convolution theorem: for any two functions f(t), g(t) of time, defined 
in the range (0, oo), their convolution is defined as: 

(f*g)= f dr f{t-r)g{r) (179) 
Jo 

Laplace transform of the (ordinary) derivative: 

n-l 

£ {/(">(*)} = s n f(s) - / 0) (0) (180) 

3=0 

where f^ n \t) is the n-th time derivative of f(t), and f(s) is the Laplace 
transform of f(t). The explicit appearance of the functions /^(0) in this 
formula makes the Laplace transform a particularly useful tool for the so- 
lution of initial value problems. We may now note that the value j = in 
the sum of the right hand side is "special": it corresponds to f^(t) = /(t), 
whereas all other terms in the sum contain true derivatives: 

n— 1 

C {/<»>(*)} = s n /(a) - s n -' /(0) - s"^ 1 / (j) (0) (181) 

3=1 

We now consider the Laplace transforms of the various objects of frac- 
tional calculus. 

a) Laplace transform of the fractional integral. 

This is easily obtained by noting that the RL integral can be viewed as a 
convolution: 

o Jf f(t) = -Sgy jT dr (t - rf- 1 /(r) = * -L f {t ) (182) 
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Hence, noting the well-known result: 

we find immediately, using the convolution theorem (179): 

£ {oJf f(t)} = 8-e f( 8 ). 

b) Laplace transform of the RL derivative. 
We write the RL derivative (155) in the form: 

D? f(t) = — g(t), m-l<P<m 

with: 

git) = JT fi fit) 

Using Eq. (180) for the m-th (ordinary) derivative, we have: 

m— 1 

C { oD? f(t)} = s m g(s) - s j g^-'HO) 

3=0 

Next, we calculate: 

Jm—j — l i ft 

g {m ~ ] ' 1] = T -T w / dr(t- r) m ^- 1 f(r) 

y dt™^- 1 V(m-(3) J v ; M ' 

Using Eq. (171), we find: 



where: 



o Jf, y. > 



From the two constraints: m — 1 < /3 < m and < j < m 
deduce: 

• when j — m — 1, we have — /? + j + 1 = m — /3 > 0, hence: (9 M 

• when j < m — 1, we have — /? + j + 1 < 0, hence: C M = -Df- 
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Note also, from Eq. (184): 



9(s)=C {oJ^ +m f(t)}=s?-™ f(s) 

We thus obtain the final form of the Laplace transform of the RL frac- 
tional derivative: 



£ {oD? f(t)} = sf> f(s)-s 

(191) 

The analogy with Eq. (181) is striking. The first term has the same form. 
In the second term the initial value of the function /(0) is replaced, for non- 
integer /3, by the initial value of the fractional integral oJf m fit), whereas 
in the terms under the sum, the initial values of the ordinary derivatives are 
replaced by the initial values of the corresponding RL fractional derivatives. 
In spite of this analogy, the RL time-derivative is not easily adapted to the 
physical applications. The Laplace transform is well known to be a valuable 
tool for the solution of initial value problems of differential equations, because 
it incorporates these initial values in Eq. (181). If we now consider initial 
value problems of fractional differential equations involving RL fractional 
derivatives with respect to time, we are faced with a serious difficulty if we 
try to use a Laplace transformation. The specification of (physical) values of 
the unknown function and of a sufficient number of its derivatives at t — 0, 
this does not allow us to determine the correponding initial valued of the 
fractional integrals and derivatives appearing in Eq. (191). This difficulty 
has been overcome by M. Caputo in [31], [32] by introducing an alternative 
definition of the fractional derivative, that is particularly useful for functions 
of a variable restricted to a finite range, such as time G (0, t). 

The Caputo fractional derivative. 

The Caputo fractional derivative is defined as follows: 



-/3+m 



f(t) 



III — Zj 



. t=0 



a — n 



oLf 3 ' 1 fit) 



t=o 



%D? fit) 



V(n-(5) 



1 d n 

(t - T y +1 - n dr™ 



f(r), n-l<(3<n. 



(192) 

Thus, the only difference with the RL derivative Eq.(156) is the commuta- 
tion of the integration and the differentiation. The Caputo derivative results 
from a composition of a fractional integral with an ordinary derivative: 



/(f) = JT fi 



d n 



(193) 
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Let us introduce a simplified notation: 

oJf-A oOt^D*, (194) 
The Caputo derivative is related to the Riemann-Liouville one as follows 

[8]: 

n-l fj-P 

D? /(f) = /(f) - E / w (0) r(j _ /g + 1) . (195) 

We give here the proof for the case < (3 < 1, which implies n = 1. The 
proof of the general case is similar, but longer. 



f(r) i r /o 



(f-r) 



/fr) 

(f-r) 

T=t 



+ dr 

T=t JO 



d_ 

dt 



(t-r)- 



dr 



1 im 



= r, m + J /TJT 1 

Hence, finally: 

DS fit) = D* /(f) - 



(t-r) J r=0 ' 7 {t-rY dr 
d 



t)P dr 



T(l-(3) 
which agrees with Eq. (195). 



t- p , 0<(3<1 



(196) 



An important property, which is obvious from Eq. (192), and also from 
(196), combined with (162), is: 

D?1 = 0, (197) 
which makes the Caputo derivative closer to the ordinary derivative. 

We now consider the composition of the Caputo derivative D% with the 
fractional integral of the same order J 13 . A first relation is immediately 
obtained from the combination of Eqs. (168) and (195): 



(198) 
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Commuting the integral and the derivative and using successively Eqs. 
(195), (169) and (152) we obtain: 



n-l 



J? /(f) = J 13 D 13 f(t) - 



t p-k 



n-l 



k=l 



r(p-k + i) 



3=0 



d j 
d? 



t=0 



and finally (with — > j + 1): 



n— 1 



^ ^ /(f) = /(*) - E { [D^- l f(t)] t=0 + /«)(o) * 



r(j + i) 

(199) 

Eqs. (198) and (199) show that the Caputo derivative is not an inverse of 
the fractional integral, either to the left or to the right. General composition 
rules of the type of Eqs. (170) - (174) can be derived for the Caputo deriva- 
tive, using the previous results; they are rather complicated, and will not be 
written down explicitly. We now consider the most important property of 
the Caputo derivative. 

Laplace transform of the Caputo derivative. 

The result is easily obtained by combining the definition (193) with Eqs. 
(184) and (181): 



ZDS fit) 



jn-p ( d n f(t) 



J t 



dt n 



n-2 



s n /(s) _ s n-l /(Q) _ J- J /("-^(O) 

and finally, setting n — k — 1 = j: 



k=0 



~ n-l 

= s 13 f(s) - s 13 " 1 /(0) - £ s ^ 3 ^ / (i) (0), n - 1< < n. 

3=1 

(200) 

This result has exactly the same form as the Laplace transform of an or- 
dinary derivative, Eq. (181), with n — > 13. Contrary to the Laplace transform 



C 



C D 13 
o 



/(f) 
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of the RL derivative (191), the Laplace transform of the Caputo derivative 
only involves the initial values of the function f(t) and of its ordinary deriva- 
tives, instead of the initial values of the fractional integral and of the RL 
fractional derivatives. In any physical differential equation, the former are 
the given input. This property makes the Caputo derivative an invaluable 
tool for solving initial value problems of fractional differential equations. 
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